classdef SimonsOriginal  <QuesNewtonBasicMethod
    properties     
    end
    methods  
        
        function obj = SimonsOriginal()
            
        end
        
          
        function [Solution obj] =Run(obj)
            %% Culot Algorithm ;
            %% A quasi-Newton algorithm for first-orderSaddle Points
            %% PDF: -> 'D:\Saddle Points Literature\Transition States\Transition
            %% State\Deterministic Methods\Single ended\A quasi-Newton algorithm for first-orderSaddle Points.pdf'

            obj.Dims = 2;
            global X
%             X = [-1; 0.0]; %% zeros(1,Dims);
            %X = [0; 0];
            
            Solution = [];
            
            %X = [-1.5; -0.5];
            %% 1) Calculate the energy and the derivatives at the initial geometry
            [obj.E, obj.Grad, obj.H]= obj.CalculateRealEnergyDerivatives(X);
            isSolution = false;            
            
            %% 2) Diagnolize the hessian matrix and choose an eigenvector
            %% vj
            if isnan(det(obj.H)) || isinf(det(obj.H)) || det(obj.H) ==0
                fprintf('\n Singular Matrix encountered during the search \n');
                return;
            end
            [obj.eigVectors eigValsM] = eig(obj.H);
            
            obj.eigVals = zeros(1,obj.Dims)
            for i = 1: obj.Dims
               obj.eigVals(i) = eigValsM(i,i);
            end    
            obj.counter = 0;            
            for K = 1: obj.MaxIterations
                %% 3) Calculate Step Vector 
                [n_pos, n_zero, n_neg] = obj.DetermineInertiaOfHessian(...
                    obj.eigVals, obj.Dims);
                obj.counter = obj.counter + 1
                Step = obj.MoveAStep();
                E1 = obj.CalculateQuadraticEnergy(Step, obj.E, obj.Grad, obj.H);    
                
                
                
                GradRMS = norm(obj.Grad);
                MaxGrad = max(obj.Grad);
                MaxStep = max(Step);
                StepRMS = norm(Step);
                
                if StepRMS == Inf
                    continue;
                end

                isSolution = obj.isATransitionState(GradRMS, MaxGrad,...
                    MaxStep, StepRMS);
                if isSolution                     
                    break;
                end

                %% 7.2.2) No:
                %% Update the geometry X, Calculate the gradient vector at the new
                %% geometry, update the hessian matrix     
                X = X + Step;
                obj.E = E1;   
                Gk = obj.Grad;
                obj.Grad = Gradient(X);
                Gk1 = obj.Grad;
                obj.H = obj.UpdateHessian(Step, Gk1, Gk, obj.H);                                  
                [obj.eigVectors eigValsM] = eig(obj.H);         
                obj.eigVals= zeros(1,obj.Dims);
                for i = 1: obj.Dims
                   obj.eigVals(i) =  eigValsM(i,i);
                end               
                
            end
            if isSolution
                Solution = X
            else
                Solution = [];
            end
            
        end
        
        function Step = MoveAStep(obj)               
            fprintf('\n Augemented QuasiNetwon Step ');            
            Step = obj.CalculateStep();            
        end
        
        function [lamba obj] = IdentifyLamba(obj)   
             d = sort(obj.eigVals);             
             factor = 0.5;
             lamba =factor*(d(1) + (d(2)/2));         %% NOTE : factor give better performance with 0.25
        end
        
        function Step = CalculateStep(obj)            
            y = obj.IdentifyLamba(); 
            Step = zeros(size(obj.Grad));
            for i = 1: obj.Dims                
                Grad_ = obj.eigVectors(:,i) * obj.Grad';               
                Grad_ = Grad_/(obj.eigVals(i)- y ) ;
                Step = Step - Grad_ * obj.eigVectors(:,i);    
            end
        end
       
    end
end
    


